#include "solve.h"
#include "matrix.h"

double* cholesky(Matrix a, double* p)       //cholesky分解法解对称正定矩阵的线性方程组
{
    a.datap[0] = sqrt(a.datap[0]);
    for (int i = 2; i < a.col+1; i++)
    {
        a.datap[i - 1] /= a.datap[0];
        a.datap[(i - 1) * a.col] = a.datap[i - 1];
    }
    for (int j = 2; j < a.row + 1; j++)
    {
        for (int k = 1; k < j; k++)
        {
            a.datap[(j - 1) * (a.col + 1)] -= a.datap[(k - 1) * a.col + j - 1]*a.datap[(k - 1) * a.col + j - 1];
        }
        a.datap[(j - 1) * (a.col + 1)] = sqrt(a.datap[(j - 1) * (a.col + 1)]);
        for (int m = j+1; m < a.col + 1; m++)
        {
            for (int k = 1; k < m-1; k++)
            {
                a.datap[(j - 1) * a.col + m-1] -= a.datap[(j - 1) * a.col + k - 1]*a.datap[(m - 1) * a.col + k - 1];
            }
            a.datap[(j - 1) * a.col + m-1] /= a.datap[(j - 1) * (a.col + 1)];
            a.datap[(m - 1) * a.col + j - 1] = a.datap[(j - 1) * a.col + m - 1];
        }
    }
    down(a, p);
    up(a, p);
    for (int w = 0; w < a.row; w++)
    {
        std::cout << p[w] <<  std::endl;
    }
}

/*int main()
{
    double b[3] = {5, 6, 5};
    double p[9] = {4, 1, 0, 1, 4, 1, 0, 1, 4};
    Matrix a(3, 3, p);
    cholesky(a, b);
}*/